#ifdef _TEST_B_
#include "Spline.h"

double B_f(double x)
{
	return 1 / (1 + 25 * x * x);
}

// f' = -50x / (1+25x^2)^2

double complete_spline(int N)
{
	INTERPC interpC;
	for (int i = 0; i < N; i++)
	{
		double x = -1 + 2.0 * i / (N - 1);
		interpC.points.push_back(x);
		interpC.interpSites[x] = {B_f(x)};
	}

	interpC.interpSites[-1].push_back(50.0 / 26 / 26);
	interpC.interpSites[1].push_back(-50.0 / 26 / 26);

	Spline<1, 3, ppForm> spline = interpolate<3, ppForm>(interpC, complete);
	//spline.plot('b');
	vector<double> err(N - 1);
	for (int i = 0; i < N - 1; i++)
	{
		double x = (interpC.points[i] + interpC.points[i + 1]) / 2;
		err[i] = spline.evaluate(x)[0] - B_f(x);
	}
	cout << spline << endl;
	return max_norm(err);
}

double notAknot_spline(int N)
{
	INTERPC interpC;
	for (int i = 0; i < N; i++)
	{
		double x = -1 + 2.0 * i / (N - 1);
		interpC.points.push_back(x);
		interpC.interpSites[x] = {B_f(x)};
	}

	Spline<1, 3, ppForm> spline = interpolate<3, ppForm>(interpC, notAknot);
	//spline.plot('b');
	vector<double> err(N - 1);
	for (int i = 0; i < N - 1; i++)
	{
		double x = (interpC.points[i] + interpC.points[i + 1]) / 2;
		err[i] = spline.evaluate(x)[0] - B_f(x);
	}
	cout << spline << endl;
	return max_norm(err);
}

void c_err(int N)
{
	cout << "re = [" << endl;
	for (int i = 10; i < N; i++)
	{
		cout << i << " " << complete_spline(i) << endl;
	}
	cout << "];" << endl;
	cout << "plot(re(:,1), re(:,2));" << endl;
}

void n_err(int N)
{
	cout << "re = [" << endl;
	for (int i = 10; i < N; i++)
	{
		cout << i << " " << notAknot_spline(i) << endl;
	}
	cout << "];" << endl;
	cout << "plot(re(:,1), re(:,2));" << endl;
}

void test_c()
{
	vector<int> num = {6, 11, 21, 41, 81};
	for (int i = 0; i < num.size(); i++)
	{
		cout << "N = " << num[i] << endl;
		cout << "Max norm: " << complete_spline(num[i]) << endl;
	}
}

void test_n()
{
	vector<int> num = {6, 11, 21, 41, 81};
	for (int i = 0; i < num.size(); i++)
	{
		cout << "N = " << num[i] << endl;
		cout << "Max norm: " << notAknot_spline(num[i]) << endl;
	}
}

int main()
{
	//cout << "x = linspace(-1, 1);" << endl;
	//cout << "plot(x, 1./(1+25*x.^2), 'r');" << endl;
	//cout << "hold on;" << endl;

	cout << "notAknot: " << endl;
	test_n();

	cout << "complete: " << endl;
	test_c();

	//c_err(200);
	//n_err(200);

	return 0;
}

#endif
